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Abstract — A Finite Element Formulation is given using the scattered field ap- 
proach with a fictitious material absorber to truncate the mesh. The formulation in- 
cludes the use of arbitrary approximation functions so that more accurate results can 
be achieved without any modification to the software. Additionally, non-polynomial 
approximation functions can be used, including complex approximation functions. 
The banded system that results is solved with an efficient sparse/banded iterative 
scheme and as a consequence, large structures can be analyzed. Results are given 
for simple cases to verify the formulation and also for large, complex geometries. 


I. Introduction 


The popularity of the Finite Element Method (FEM) in Electromagnetic Scattering is largely 
due to the fact that the method results in a banded system with an 0(N) memory demand. This 
is in sharp contrast to the Moment Method which has an 0(A 2 ) memory demand. For open 
region scattering problems however, the Finite Element region must be artificially truncated in 
order to bound the problem. This has been accomplished with approximate methods such as 
the Uni-Moment Method [1] and the Absorbing Boundary Condition [6] and more successfully 
with the rigorous Boundary Integral Technique [3]. Unfortunately, the solution of the Boundary 
Integral in general results in a dense system thus destroying the 0(N) character of the FEM. 
Also, the boundary integral technique breaks down at resonant frequencies and these resonant 
frequencies become more and more closely spaced as the size of the boundary is increased. 

More recently, the technique of using a fictitious material absorber to confine the scattered 
field has been suggested in [4]. Since fictitious materials having a negative real part of \i and 
f can be used, an absorber with a very low reflection coefficient even at near grazing incidence 
can be realized. Also, this technique requires the unknown quantity to be the scattered field, 
not the total field as is usually the case. The benefits of this are two-fold. First, the scattered 
field does not vary as rapidly as the total field and second, the scattered field approach requires 
that the incident field be evaluated in the interior of the mesh, thus reducing the phase error 
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Figure 1: Region of solution 


normally associated with the total field approach. That is, the total field approach requires that 
the incident field be evaluated only at the fictitious boundary and then propagate through the 
mesh to the scatterer thus introducing phase error, especially notable when the mesh is large. It 
has been confirmed in [5] that the scattered field approach results in a more accurate solution. 

It is known that there is also an overall inherent error in the FEM that does not exist in the 
Moment Method [7]. This error is a result of the bandedness of the system and is introduced 
as the solution propagates through the mesh. In fact, experience has shown that as the mesh 
area increases, the element size must decrease. The only solution to this problem is to use higher 
order approximation functions, which usually involves a complete re-formulation and re-design 
of the software. 

In this paper a number of advances in FEM calculations are combined to give a solution 
to the 2-Dimensional Helmholtz equation for electrically large scatterers. The scattered field 
approach is used and a fictitious material absorber is used to truncate the mesh. Unique to this 
formulation is the use of arbitrary approximation functions that make it possible to increase the 
accuracy of the solution without modifying the software. 


II. Formulation 


Consider the region shown in Figure 1. 

Given that there is some incident field, we wish to solve the wave equation in the region Q 
and find the fields scattered by the dielectric/magnetic/conducting target. 
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The source free wave equation in two dimensions can be written 


V • + klv<j> = 0 


where 


<t> = E : 

U = 1//J r 
V = 1/fr 


for E polarization ( z directed electric field) 
and 


<t>=H; 

u = l/e r 

V = 1 /fi r 

for //polarization (z directed magnetic field) 


( 1 ) 


( 2 ) 

(3) 

(4) 


(5) 

(6) 
(7) 


Now let T be some compactly supported testing function used to enforce (1) over the entire 
region . 



T[V ■ (uV</>) 


+ klv<t>] ds = 0 


( 8 ) 


Using the identity 


T[V • (t/Vtf)] = -u(V0 ■ Vr) + V • (TuX?<j>) 


(9) 


gives 



[u(VT • V<?!>) 



V • (TuV<t>)ds — 0 


( 10 ) 


and applying the divergence theorem to the second term in the above brings us to the weak form 
of eqn . ( 1 ) . 


J J [u(VT • V<f>) - k\vT4>] ds — <j> uT di - j> uT ^ d£ = 0 




(ID 


To use eqn. (11) in an FEM implementation would be considered a total field approach since 
the unknowns would be the total field. To get the analog of (11) with the scattered field as the 
unknown we let 


<t> = <b s 4- <f>' 


( 12 ) 
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and (11) becomes 


Jj[u<?T • V** ) - klvTtf] ds- £ T(u?£- + a ?£)dt 
£ T(u^ + u^£-)d(= - J £[u(VT -W) - kZiTd>']ds 


where il and v are the material constants of the scatterer in a free space environment while 
u and v include the fictitious absorber. In this way the incident field does not interact with the 
fictitious absorber. 

In either case, eqns. (11) and (13) require that the normal derivative of the unknown (total 
field and scattered field respectively) be specified on the enclosing boundaries r a and r c . These 
normal derivatives are arbitrary, and once given, a solution in the interior is uniquely defined. 

For E polarization we will put a perfect magnetic conductor at T a and for H polarization we 
will put a perfect electric conductor at T a thus the second term on the LHS of (13) is 




Note that this conductor is seen by the scattered field only, not the total field. 


The boundary condition at T c is 


<9n <9n 


for H polarization, and thus the third term on the LHS of (13) is zero. 

For E polarization the boundary condition at T c is 

^ + ^ = 0 (16) 

and since the solution is specified on this contour, the normal derivatives will uncouple from the 
rest of the equation once we have discritized (13) and created a linear system, thus the third 
term on the LIIS of (13) drops out. 

We have that eqn. (13) becomes 

J jf[l//i,(vr -W) - lcl< r T<t>‘]ds - £ 1 fr rT ?£dt= (17) 

- JJll/WT ■ W) - J fcg W] ds 

for E polarization and 

/ / n [l/c r (VT • W‘) - *S/i r 2V] ds - £ 1 /I t Tj£- d£ = (18) 

- Jj{l/lr(VT ■ W) - kliir T4?) ds 
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for polarization. 

The region Q now is divided into arbitrary, polygonal finite elements. Each element is defined 
by its corner nodes and may contain other nodes located anywhere in the element. The solution 
in the e’th element is approximated by 


N e 

<t> e (x,y) = ^S'(x, (19) 

n = l 

where N e is the total number of nodes in the e’th element, <f> e n is the value of <fi at the n’th node 
and S* (x, y) is the n’th shape function for the e’th element understood to be zero outside of the 
e’th element. These shape functions will be defined in the next section. To solve (17) and (18) 
using Galerkin’s method, we let the testing function T be the same as the basis function in (19). 


Substituting (19) into (17) and (18) gives the Element Equations 

/ f [1/tfVSf ■ V Y, W ~ ^ £ Sj 4>)) ds = 

J •'O. j - 1 j = 1 

+ jf 1/#S? J^ d( ~ JJ Q U/#(VS? • Vd»‘) - kge‘S- 4>') ds 


for E polarization and 


i = 1 , 2 , 3 . . . A r 


I ds — 


( / [1 A*V5? • V J2 W - k lSr s t E s 'y W 1 

J j = 1 > = l 

+ jf i/V r Sf^di- [1/$(VS? .VW-klfcSfWds t = l,2,3..JV 


for H polarization. 


( 20 ) 


( 21 ) 


Again, jl* and c e r are the material constants of the scatterer in a free space environment while 
and e* include the fictitious absorber. 


Equations (20) and (21) define an A rc by N e linear system for a single element. Assembling 
all of the element systems together and converting to global node numbers gives the global finite 
element system 

[*]{*'} = {**} ( 22 ) 


For £ polarization only, the boundary condition (16) must be enforced by setting <f> 3 = — <f> % 
for nodes on the conductor. 


III. Arbitrary Approximation Functions 
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In eqn. (19) the expansion of the solution was given as a sum of A f shape functions weighted 
by the nodal values of the solution. To determine these shape functions let 


N m 


<t> e {x,y) = ^a n / n (x,y) 

n = l 


(23) 


where the / n (x,y)’s are linearly independent functions known as approximation functions. 

Forcing <j> c (x,y) to be 4>*(xi , y t ) at each nodal point gives the system 

t^]K} = m (24) 

where [T£] is an N by N matrix whose i.j’th element is given by fj (x, , y t ) and x^y, are 
the coordinates of the Vth node. For ease of programming, these approximation functions are 
often chosen to be linear, ie. f\ = l,/ 2 = x ,/3 = y. However, since it is known that higher 
order approximation functions result in better conditioned systems and are more accurate 1 , it is 
desirable to use approximation functions of arbitrary order or even of a non-standard type such 
as complex transcendental functions which are better suited for Electromagnetics. 

The shape functions are then given by 

N m 

Sn(x>v) = '52fi(*,y)[F' n ]- 1 (25) 

1=1 

Notice that the first derivatives of the shape functions are also required and this must be taken 
into account when choosing approximation functions since from (23), derivatives of the shape 
function transfer directly to the approximation functions. 

Note also that the inverse of an N € by N e matrix is required to calculate the shape function 
for each element. Since N e will not be very large, this does not introduce any great expense. 

When the arbitrary shape functions are used to generate the element matrices an efficient 
numerical integration scheme must be used. Since the elements are arbitrary polygons, an in- 
tegration algorithm based on the polygon fill algorithm from computer graphics w T as developed. 
This integration algorithm takes as input: the number of points defining the polygon, the coordi- 
nates of these points, the number of sub-regions to split the area into, the order of the Gaussian 
quadrature to use, and the integrand function. With some experimentation, it was found that 
increasing the accuracy of the integrations had little effect on the far field scattering patterns 
once the integrations were approximately five percent accurate (for linear shape functions.) 


IV. Results 


The bistatic RCS pattern for a perfectly conducting circular cylinder of radius A/2 is shown in 
figure 2. Linear shape functions were used and as can be seen from figure 3 a discretization 
of A/20 was required. The fictitious absorber was placed A/2 away from the cylinder to avoid 
interactions with the scatterer. 
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Figure 2: Bistatic RCS of conducting cylinder 



Figure 3: Mesh for r = A/2 conducting cylinder 
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Figure 4: Bistatic RCS of Airfoil 


A more useful result is shown in figure 4 for a perfectly conducting airfoil at low frequency. 
The wing is approximately if A long and as can be seen from figure 5 a discretization of A/20 
was again required but the fictitious absorber needed to be placed a full wavelength away from 
the target in order to get accurate results for H polarization. The reason for the discrepancy in 
the E polarization results is that the artificial absorber actually needs to be placed farther away 
from the scatterer. The following discussion addresses this issue. 

For quadratic shape functions, which are obtained by setting /i = l,/ 2 = x,/ 3 = y, / 4 = 

/s = x 1 2 , and / 6 = y 2 , the mesh is not as dense as can be seen from figure 6. However since 
there are twice as many nodes per element, the net effect is that there are about half as many 
elements but the same number of nodes. However the increase in accuracy makes it possible to 
analyze larger and more detailed scatterers. 

For cubic shape functions, which are obtained by setting, for example, f\ = l,/ 2 = x,/ 3 = 
y,/ 4 = xy, / 5 = x 2 ,/ 6 — y 2 ,/7 = x 3 ,/ g = y 3 ,/^ = a; 3 */ 3 , the mesh is even less dense as can be 
seen from figure 7. Again, there are approximately the same number of nodes but the number 
of elements is cut in half. Actually, since the accuracy is increased exponentially, the mesh size 
could be increased further thus actually reducing the number of unknowns. This would seem 
to imply that higher order approximation functions keep getting more efficient. Whether or not 
this is the case remains to be seen. 

Another choice of approximation function could be the Hankel function of order 0. This would 
make it possible to discretize using larger elements far away from the scatterer where the wave 
becomes somewhat cylinderical. There would still be other approximation functions to allow' 

1 Polynomial approximation functions of Order p result in a solution which is h p accurate where h is the mesh 

size. 
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Figure 5: Mesh for Airfoil : Linear 



Figure 6: Mesh for Airfoil : Parabolic 
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Figure 7: Mesh for Airfoil : Cubic 


some modification to the cylinderical wave. Being able to efficiently discretize large regions is 
desirable since it is clear from the results for the airfoil that the artificial absorber must be placed 
far away from the scatterer. Actually, one could consider the transcendental Hankel function to 
be of exact accuracy limited only by numerical considerations. This is a feature that polynomial 
basis functions do not possess since they do not themselves solve the wave equation. 


V. Conclusions 


Although much more testing must be done with both the placement of the artificial absorber 
and the use of higher order and complex approximation functions, it is clear that this type of 
formulation holds promise for the accurate analysis of large scattered. While many formulations 
are specific to a certain problem and do not allow for growth and modification, the above formu- 
lation and the software developed provide a flexible tool for improving the current state of FEM 
for the analysis of large electromagnetic scattering problems. 
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